R jako GIS
Podpůrný text pro předmět: GIS pro biologické aplikace
Autor: Matěj Man
Aktualizace: 08. 11. 2021
Odkaz na R skript: Kopírovat a vložit do R
Odkaz na data: Stáhnout data zip Odkaz na github repozirář: Repozitář
Další zdroje:
- Online kniha Geocomputation with R
- Online kniha Spatial Data Science with R
- Velice pěkný návod pro R knihovnu sf
- The Visual Raster Cheat Sheet Rpubs dokument
Knihovny
# list.of.packages <- c("sf","raster","mapview","whitebox","randomcoloR","leaflet","Rcpp")
# install.packages(list.of.packages)
# install.packages("devtools")
# install.packages('raster', repos='https://rspatial.r-universe.dev')
# library(devtools)
# install_github("r-spatial/sf")
library(sf)
library(raster)
library(mapview)
library(randomcoloR)
library(leaflet)Načítání vektorových GIS dat do R
## Načítání vektorových dat shp
# nastavte kde máte u sebe na PC data
# pozor musíte zdvojit nebo otočit lomítka
cesta<-"d:/Git/gisproba/data/"
# zkonstuuje cestu kde leží vrstva Brdy
data.path<-paste0(cesta,"CHKO_Brdy.shp")
# načte .shp do R
brdy<-st_read(data.path,stringsAsFactors = F,quiet = T)
# obdobně načteme třeba hranici ČR
data.path<-paste0(cesta,"hrcr1_wgs.shp")
hrcr<-st_read(data.path,stringsAsFactors = F,quiet = T)
# prostý obrázek, bez interaktivity
plot(st_geometry(hrcr))
# parametr add přidá další vrstvu do existujícího obrázku
plot(st_geometry(brdy),col="red",add=T) 
# Nastaví pracovní adresář working dorectory (WD)
# dále nemusíme data z WD volat celou cestou, stačí názvy
setwd(cesta)
## Načítání vektorových dat z tabulky
# načíst prostou tabulku
chmu<-read.table("staniceCHMUtablecoma.csv",header = T, sep=",")
# převést tabulky na prostorová data
chmu_sf<-st_as_sf(chmu, coords = c("Xcoo", "Ycoo"),crs = 4326)
## vizualizovat prostorová data interaktivně
# pokročilejší balíček leaflet
# definujeme paletu o 17 barvách podle typu stanice
distCol<- colorFactor(distinctColorPalette(17), chmu_sf$Station.ty)
# definoujeme okno s mapou
leaflet(chmu_sf) %>%
addTiles() %>%
addCircleMarkers(popup=chmu_sf$Name,color = ~distCol(Station.ty),
stroke = FALSE, fillOpacity = 0.8,radius=4) %>%
addLegend(pal = distCol, values = ~Station.ty)Projekce a crs
st_crs(brdy) # jaký crs má vrstva brdy ?## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]st_crs(hrcr) # jaký crs má vrstva hrnice čr ?## Coordinate Reference System:
## User input: WGS 84 + Unknown VCS from ArcInfo Workstation
## wkt:
## COMPOUNDCRS["WGS 84 + Unknown VCS from ArcInfo Workstation",
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["Degree",0.0174532925199433]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["Degree",0.0174532925199433]]],
## VERTCRS["Unknown VCS from ArcInfo Workstation",
## VDATUM["Unknown"],
## CS[vertical,1],
## AXIS["gravity-related height (H)",up,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]]## vektory
# transformace podle EPSG
brdy_32633<-st_transform(brdy,32633)
st_crs(brdy_32633) # jaký crs má vrstva brdy_32633 ?## Coordinate Reference System:
## User input: EPSG:32633
## wkt:
## PROJCRS["WGS 84 / UTM zone 33N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 33N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",15,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Between 12°E and 18°E, northern hemisphere between equator and 84°N, onshore and offshore. Austria. Bosnia and Herzegovina. Cameroon. Central African Republic. Chad. Congo. Croatia. Czechia. Democratic Republic of the Congo (Zaire). Gabon. Germany. Hungary. Italy. Libya. Malta. Niger. Nigeria. Norway. Poland. San Marino. Slovakia. Slovenia. Svalbard. Sweden. Vatican City State."],
## BBOX[0,12,84,18]],
## ID["EPSG",32633]]# transformace podle proj4 string
brdy_5514<-st_transform(brdy, "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56")
st_crs(brdy_5514) # jaký crs má vrstva brdy_5514 ?## Coordinate Reference System:
## User input: +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56
## wkt:
## BOUNDCRS[
## SOURCECRS[
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["Unknown based on Bessel 1841 ellipsoid",
## ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]],
## CONVERSION["unknown",
## METHOD["Krovak (North Orientated)",
## ID["EPSG",1041]],
## PARAMETER["Latitude of projection centre",49.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8811]],
## PARAMETER["Longitude of origin",24.8333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8833]],
## PARAMETER["Co-latitude of cone axis",30.2881397222222,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",1036]],
## PARAMETER["Latitude of pseudo standard parallel",78.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8818]],
## PARAMETER["Scale factor on pseudo standard parallel",0.9999,
## SCALEUNIT["unity",1],
## ID["EPSG",8819]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Position Vector transformation (geog2D domain)",
## ID["EPSG",9606]],
## PARAMETER["X-axis translation",570.8,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",85.7,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",462.8,
## ID["EPSG",8607]],
## PARAMETER["X-axis rotation",4.998,
## ID["EPSG",8608]],
## PARAMETER["Y-axis rotation",1.587,
## ID["EPSG",8609]],
## PARAMETER["Z-axis rotation",5.261,
## ID["EPSG",8610]],
## PARAMETER["Scale difference",1.00000356,
## ID["EPSG",8611]]]]# transformace podle jiné existující vrstvy
brdy_krovak<-st_transform(brdy, st_crs(brdy_5514))
st_crs(brdy_krovak) # jaký crs má vrstva brdy_5514 ?## Coordinate Reference System:
## User input: +proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56
## wkt:
## BOUNDCRS[
## SOURCECRS[
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["Unknown based on Bessel 1841 ellipsoid",
## ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]],
## CONVERSION["unknown",
## METHOD["Krovak (North Orientated)",
## ID["EPSG",1041]],
## PARAMETER["Latitude of projection centre",49.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8811]],
## PARAMETER["Longitude of origin",24.8333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8833]],
## PARAMETER["Co-latitude of cone axis",30.2881397222222,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",1036]],
## PARAMETER["Latitude of pseudo standard parallel",78.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8818]],
## PARAMETER["Scale factor on pseudo standard parallel",0.9999,
## SCALEUNIT["unity",1],
## ID["EPSG",8819]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]],
## TARGETCRS[
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]],
## ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
## METHOD["Position Vector transformation (geog2D domain)",
## ID["EPSG",9606]],
## PARAMETER["X-axis translation",570.8,
## ID["EPSG",8605]],
## PARAMETER["Y-axis translation",85.7,
## ID["EPSG",8606]],
## PARAMETER["Z-axis translation",462.8,
## ID["EPSG",8607]],
## PARAMETER["X-axis rotation",4.998,
## ID["EPSG",8608]],
## PARAMETER["Y-axis rotation",1.587,
## ID["EPSG",8609]],
## PARAMETER["Z-axis rotation",5.261,
## ID["EPSG",8610]],
## PARAMETER["Scale difference",1.00000356,
## ID["EPSG",8611]]]]# kde jsou ty brdy? no nevidím je protože to má jiný CRS
plot(st_geometry(hrcr),main="Kde jsou ty Brdy?")
plot(st_geometry(brdy_32633),add=T) 
# musíme tedy transformovat jedno nebo druhé
hrcr_32633<-st_transform(hrcr, 32633)
plot(st_geometry(hrcr_32633),main="No jo, už máme správný CRS")
plot(st_geometry(brdy_32633),add=T,col="red")
Načítání rastrových GIS dat do R
# načítání jdnoduchou funkcí "raster"
DEM<-raster(paste0(cesta,"DEM_Jested_100m.tif"))
# kontrolí obrázek
# data tu jsou
plot(DEM)
# ale leží na správném místě?
# ukáže crs vrstvy
DEM@crs## CRS arguments:
## +proj=krovak +axis=swu +lat_0=49.5 +lon_0=24.8333333333333 +alpha=0
## +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs# zkušeným okem, nebo podle errorů při další alanýze poznáte,
# že nemá dobře definovaný crs, je to křovák, ale šptně definovaný.
# tak mu řekneme, jaká je korektní definice křováka
# http://freegis.fsv.cvut.cz/gwiki/S-JTSK
DEM@crs<-crs("+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +alpha=30.28813972222222 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +pm=greenwich +units=m +no_defs +towgs84=570.8,85.7,462.8,4.998,1.587,5.261,3.56")
# i rastr je možné vizualizovat interaktivně
# definujeme barvy
pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(DEM),
na.color = "transparent")
leaflet(DEM) %>%
addTiles() %>%
addRasterImage(DEM,colors = pal, opacity = 0.8) %>%
addLegend(pal = pal, values = values(DEM),
title = "Elevation [m]")Příklad vektorové analýzy
# plocha polygonů
st_area(brdy)## 343898692 [m^2]# délka linií (obvod polygonů)
st_length(brdy)## 0 [m]# buffer 5 km
brdy_5000<-st_buffer(brdy_32633,5000)
plot(st_geometry(brdy_5000), main="buffer 5 km")
plot(st_geometry(brdy_32633),add=T,col="red")
# centroidy
brdy_cen<-st_centroid(brdy_32633) ## Warning in st_centroid.sf(brdy_32633): st_centroid assumes attributes are
## constant over geometries of xplot(st_geometry(brdy_32633), main="centroid",graticule=T,axes=T)
plot(st_geometry(brdy_cen),add=T,col="blue",pch=19)
# průnik
chmu_brdy<-st_intersection(brdy,chmu_sf)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesplot(st_geometry(brdy), main="intersection - chmu stanice v brdech",graticule=T,axes=T)
plot(st_geometry(chmu_brdy),add=T,col="red",pch=19)
Hromadné zpracování dat
Aneb co by nám v QGIS trvalo klikat hodiny můžeme v R naprogramovat během minut
# načteme všechny cesty k shp vrstvám v daném adresáři
parky.files<-list.files(path=cesta,pattern="*.shp$",full.names = T)
parky.files<-parky.files[-1]
## když chceme všechny soubor do jednoho objektu
# načte první shpfile
parky.init<-st_read(parky.files[1],quiet = TRUE)
# připojí všechny další shp file do jendoho objektu
for (i in 2:length(parky.files)) {
parky.next<-st_read(parky.files[i],quiet = TRUE)
parky.init<-rbind(parky.init,parky.next)
}
head(parky.init)## Simple feature collection with 6 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 13.6537 ymin: 48.81219 xmax: 18.74347 ymax: 50.67169
## Geodetic CRS: WGS 84
## OBJECTID geometry
## 1 3900 POLYGON ((18.53572 49.66022...
## 2 3894 POLYGON ((17.88455 49.15734...
## 3 3873 POLYGON ((14.82518 49.66756...
## 4 3877 POLYGON ((14.19356 49.00166...
## 5 3904 POLYGON ((13.93738 49.81539...
## 6 3888 POLYGON ((16.10584 50.662, ...# připojit metadat GIS join
meta.path<-paste0(cesta,"metadata_parky.csv")
metadata<-read.table(meta.path,sep=";",header=T,stringsAsFactors = T)
parky<-merge(parky.init,metadata)
head(parky)## Simple feature collection with 6 features and 10 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 13.61426 ymin: 48.80132 xmax: 15.00087 ymax: 50.67842
## Geodetic CRS: WGS 84
## OBJECTID KOD KAT NAZEV ROZL ZMENA_G ZMENA_T
## 1 3873 21 CHKO Blaník 4029.195 20030101 20170901
## 2 3874 22 CHKO Český kras 13225.701 20120824 20170901
## 3 3875 23 CHKO Kokořínsko - Máchův kraj 41037.143 20150223 20170901
## 4 3876 24 CHKO Křivoklátsko 62497.146 20131107 20170901
## 5 3877 31 CHKO Blanský les 21968.998 20091120 20170901
## 6 3878 32 CHKO Třeboňsko 68744.620 20091218 20170901
## PREKRYV SHAPEAREA SHAPELEN geometry
## 1 0 40291954 31958.58 POLYGON ((14.82518 49.66756...
## 2 0 132257006 83895.23 POLYGON ((14.3244 50.00705,...
## 3 0 410371433 198507.51 MULTIPOLYGON (((14.4139 50....
## 4 0 624971460 138706.87 POLYGON ((13.81561 50.14976...
## 5 0 219689976 79790.35 POLYGON ((14.19356 49.00166...
## 6 0 687446196 143605.35 POLYGON ((14.71406 49.18412...# kontrolní obrázek
plot(st_geometry(hrcr))
plot(st_geometry(parky),add=T,col="green")
# interaktivně
distCol<- colorFactor(distinctColorPalette(32), parky.init$OBJECTID)
leaflet(parky) %>%
addTiles() %>%
addPolygons(color = ~distCol(OBJECTID),
stroke = FALSE, fillOpacity = 0.8,
label = paste0(parky$KAT,"_",parky$NAZEV))Rastrové analýzy
Velice dobrá knihovna pro rastrové analýzy v R je whitebox, což je implementace jinak v Python/Rust napsaných geo algoritmů.
## instalace balíčku whitebox pro rastrové analýzy
# install.packages("whitebox", repos="http://R-Forge.R-project.org") # instalace
# whitebox::wbt_init() # inicializace
library(whitebox)
# pozor na diakritiku a mezery v cestě. Pro WBT nesmí být
# cesta k rastrovým datům
dem <-("d:/Git/gisproba/data/DEM_Jested_100m.tif")
## stínovaný reliéf
# cesta k budoucímu výsledku
output<- ("c:/Users/Public/Documents/Hillshade_100m.tif")
# spustí algoitmus pro výpočet stinovaného reliéfu
wbt_hillshade(dem, output, azimuth = 315, altitude = 30,
zfactor = 1,verbose_mode = FALSE)## [1] "hillshade - Elapsed Time (excluding I/O): 0.1s"# načtu výsledek jako rastr
hill<-raster(output)
# načtu původní dem jako rastr
elev<-raster(dem)
# vytisknu oba obrázky pro kontrolu
plot(hill, col=grey(0:100/100), legend=FALSE,main="stínovaný reliéf")
plot(elev, col=rainbow(25, alpha=0.35), add=TRUE)
## Sklon svahu
output<- ("c:/Users/Public/Documents/slope_100m.tif")
wbt_slope(dem,output,zfactor = 1, verbose_mode = FALSE)## [1] "slope - Elapsed Time (excluding I/O): 0.1s"slp<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="Sklon svahu")
plot(slp, col=heat.colors(25, alpha=0.2), add=TRUE)
## The terrain ruggedness index (TRI)
output<- ("c:/Users/Public/Documents/TRI_100m.tif")
wbt_ruggedness_index(dem,output,zfactor = 1, verbose_mode = FALSE)## [1] "ruggedness_index - Elapsed Time (excluding I/O): 0.1s"tri<-raster(output)
plot(hill, col=grey(0:100/100), legend=FALSE,main="terrain ruggedness index")
plot(tri, col=cm.colors(25, alpha=0.3), add=TRUE)